%% Load the data - Select the File: 'SwimmingSpeedVsBulkViscosity'
close all; clear all; clc;
clearvars;
total_directory=uigetdir(); %pick the folder you want to work in 
total_directory = dir(total_directory);
dirFlags = [total_directory.isdir];
total_directory = total_directory(dirFlags);
total_directory = total_directory(~ismember({total_directory.name},{'.','..'}));
%% Filter Tracks
mag = 10; fps = 30;Camera_um_pix = 6.5;
TracksSel = [];YwallMax_s = [];YwallMin_s = [];myparameters = [];
myparameters.lengthdev = 4;
h = waitbar(0,'Filtering');
for i = 1:length(total_directory)
    cd([total_directory(i).folder filesep total_directory(i).name]);load('ParticleTracks');
    [tracks_c] = Camera2MicronSecond(tracks,mag,fps,Camera_um_pix);
    [tracks_filt, trackstats] = FilterTracks(tracks_c, myparameters);
    MeanSpeed(i) = mean(vertcat(tracks_filt(:).Speed_um));Std(i) = std(vertcat(tracks_filt(:).Speed_um));
    TracksSel{i} = tracks_filt;
    len(i) = length(tracks_filt);
    clear tracks tracks_c tracks_filt
    waitbar(i / length(total_directory),h)    
end
close(h);'done'
OverallMean = [(mean([MeanSpeed(1) MeanSpeed(2)])) (mean([MeanSpeed(3) MeanSpeed(4) MeanSpeed(5)]))...
    (mean([MeanSpeed(6) MeanSpeed(7) MeanSpeed(8)])) (mean([MeanSpeed(9) MeanSpeed(10) MeanSpeed(11)]))...
    (mean([MeanSpeed(12) MeanSpeed(13) MeanSpeed(14)])) (mean([MeanSpeed(15) MeanSpeed(16)]))];
MeanStandardDeviation = [mean([Std(1) Std(2)]) mean([Std(3) Std(4) Std(5)])...
    mean([Std(6) Std(7) Std(8)]) mean([Std(9) Std(10) Std(11)])...
    mean([Std(12) Std(13) Std(14)]) mean([Std(15) Std(16)])];
%% Viscosity Calculation
PerPEO = [0 0 0.125 0.125 0.125 0.25 0.25 0.25 0.375 0.375 0.375 0.5 0.5 0.5 0.82 0.82];
Percentage_ref = [0 0.25 0.5 1.0];
Viscosity_ref = [1 2 5 12];
p = polyfit(Percentage_ref, Viscosity_ref,2);x_eta = linspace(0,1,100);
p_fit = p(3)+p(2).*x_eta+p(1).*x_eta.^2;

con = [0 0.125 0.25 0.375 0.5 0.82];
Viscosity_calc = (p(3)+p(2).*con + p(1).*con.^2)
Viscosity = [0.9959 Viscosity_calc(2) 2.2341 Viscosity_calc(4) 4.4814 8.9618]
%% Fit a line to the Swimming Speed vs Bulk Viscosity
close all; clc;
alw = 0.5;  % Border line width
afs = 10;   % Font Size
LW = 0.5;   % Error Bar Line Width
clw = 1;    % Plot Line Width
Cap = 0;
markersize = 3.5;
mypapersize = [1.53 1]; % [width height]
LineColor = [0 0 0];

x_eta = linspace(0.82,10,100);
Const = polyfit(log(Viscosity),log(OverallMean), 1);
m = Const(1);k = Const(2);
V_fit = x_eta.^m.*exp(k);

triang_x = [2.5, 3.4];  % location for slope triangle
xx = linspace(triang_x(1),triang_x(2),10)
yy = 220*xx.^-1
triang_y = interp1(xx, yy, triang_x);

SpeedvsViscosity_fit = figure('color',[1 1 1]);
    plot(x_eta,V_fit,'color',LineColor,'LineWidth',LW);
    hold on;
    errorbar(Viscosity,OverallMean,MeanStandardDeviation,'o','CapSize',Cap,...
        'Color',LineColor,'MarkerEdgeColor',LineColor,'LineWidth',clw,'MarkerSize',...
        markersize, 'MarkerFaceColor','w');
    plot(triang_x([1,2,2,1,2]), triang_y([1,1,2,1,2]),'color',LineColor,'LineWidth',LW)
    hold off;
    set(gca,'YScale','log','XScale','log');
    xlim([0.9 10]);ylim([0 152]);
    box on;
    xticks([0,1,5,9]);
    xticklabels({'','','','',''});
    yticks([0,20,50,100]);
    yticklabels({'','','',''});
%
    set(gca,'XColor','k');set(gca,'YColor','k');set(gca,'LineWidth',alw)
    set(gcf, 'PaperUnits', 'inches', 'papersize', mypapersize, 'PaperPosition', [0 0 mypapersize]);
    set(SpeedvsViscosity_fit,'Units','inches','Position',[0 0 mypapersize]);
%
    ax = gca;outerpos = ax.OuterPosition; ti = ax.TightInset;
    left = outerpos(1);bottom = outerpos(2);ax_width = outerpos(3);ax_height = outerpos(4);
    ax.Position = [left+0.009 bottom+0.015 ax_width-0.025 ax_height-0.03];
%% Save Figure
 figdir = '.\Figure1';
 print(SpeedvsViscosity_fit, '-painters','-dpdf', '-r1200', [figdir,'V6_SD.pdf']); 
%% Save Data in Excel
dir_Excel = '.\Source data\MainText';
writematrix(Viscosity',[dir_Excel,'Figure1.xlsx'],'Sheet',6,'Range','A3')
writematrix(OverallMean',[dir_Excel,'Figure1.xlsx'],'Sheet',6,'Range','B3')
writematrix(MeanStandardDeviation',[dir_Excel,'Figure1.xlsx'],'Sheet',6,'Range','C3')
writematrix(x_eta',[dir_Excel,'Figure1.xlsx'],'Sheet',6,'Range','D3')
writematrix(V_fit',[dir_Excel,'Figure1.xlsx'],'Sheet',6,'Range','E3')
%% FUNCTION: Convert pix/frame to um/s
function [c_tracks] = Camera2MicronSecond(tracks,mag,fps,Camera_micron_pix)
% Camera2MicronSecond: Converts camera units (frames and pix) to physical
% units (microns and seconds). Cameras have different pixel dimensions. 
% For example, the Zyla has 6.5 microns per pixel
c_tracks = tracks;

for i=1:length(c_tracks)
    c_tracks(i).t_sec = c_tracks(i).t/fps;
    c_tracks(i).x_um = c_tracks(i).x.*Camera_micron_pix/mag;
    c_tracks(i).y_um = c_tracks(i).y.*Camera_micron_pix/mag;
    c_tracks(i).u_um = c_tracks(i).u.*Camera_micron_pix.*fps./mag;
    c_tracks(i).v_um = c_tracks(i).v.*Camera_micron_pix.*fps./mag;    
    % Calculate the Speed of the Trajectory (magnitude of u-velocity and v-velocity)    
    s = sqrt(c_tracks(i).v.^2+c_tracks(i).u.^2);
    c_tracks(i).Speed = s;
    
    s_um = sqrt(c_tracks(i).v_um.^2+c_tracks(i).u_um.^2);
    c_tracks(i).Speed_um = s_um;   
end
end
%% FUNCTION: filter tracks
function [f_tracks, trackstats] = FilterTracks(tracks, myparameters)

if ~isfield(myparameters,'minmaxlengthtime')||isempty(myparameters.minmaxlengthtime); myoptions.minmaxlengthtime=0; else; minlengthtime = myparameters.minmaxlengthtime(1);maxlengthtime=myparameters.minmaxlengthtime(2); myoptions.minmaxlengthtime=1; end
if ~isfield(myparameters,'minmedianspeed')||isempty(myparameters.minmedianspeed);  myoptions.minmedianspeed=0; else; minmedianspeed = myparameters.minmedianspeed; myoptions.minmedianspeed=1; end
if ~isfield(myparameters,'minRG') || isempty(myparameters.minRG); myoptions.minRG = 0; else; minRG = myparameters.minRG; myoptions.minRG = 1; end
if ~isfield(myparameters,'lengthdev')||isempty(myparameters.lengthdev); myoptions.lengthdev=0; else; lengthdev=myparameters.lengthdev; myoptions.lengthdev=1; end
if ~isfield(myparameters,'minmaxcontourlength')||isempty(myparameters.minmaxcontourlength);myoptions.minmaxcontourlength=0; else; mincontourlength=myparameters.minmaxcontourlength(1); maxcontourlength=myparameters.minmaxcontourlength(2); myoptions.minmaxcontourlength=1; end

f_tracks=tracks;

if myoptions.minmaxlengthtime==1
    ind_minmaxlengthtime=zeros(1,length(f_tracks));tracklength=zeros(1,length(f_tracks));
    for k=1:length(f_tracks)
        tracklength(k) = length(f_tracks(k).x);
        if tracklength(k) < maxlengthtime && tracklength(k) > minlengthtime
            ind_minmaxlengthtime(1,k)=1;
        end
    end
    trackstats.tracklengthtime=tracklength;
else
    ind_minmaxlengthtime=ones(1,length(f_tracks));
end

if myoptions.minmedianspeed==1
    med=zeros(1,length(f_tracks));ind_minmedianspeed=zeros(1,length(f_tracks));
    for k=1:length(f_tracks)
        u_r = sqrt(tracks(k).u.^2 + tracks(k).v.^2);
        med(k) = median(u_r);
        if med(k) > minmedianspeed
            ind_minmedianspeed(1,k)=1;
        end
    end
    trackstats.medianspeed=med;
else
    ind_minmedianspeed=ones(1,length(f_tracks));
end

if myoptions.minRG==1
    ind_minRG=zeros(1,length(f_tracks));RG=zeros(1,length(f_tracks));
    for k=1:length(f_tracks)
        RG(1,k)= sqrt((max(f_tracks(k).x)-min(f_tracks(k).x))^2+(max(f_tracks(k).y)-min(f_tracks(k).y))^2);
        if RG(1,k)>minRG
            ind_minRG(1,k)=1;
        end
    end
    trackstats.RG=RG;
else
    ind_minRG=ones(1,length(f_tracks));
    
end

if myoptions.lengthdev==1
    ind_lengthdev=zeros(1,length(f_tracks));RG=zeros(1,length(f_tracks));contourlength=zeros(1,length(f_tracks));dev=zeros(1,length(f_tracks));
    s_tracks=EqualSpace(f_tracks,1);
    for k=1:length(f_tracks)
        RG(1,k)= sqrt((max(f_tracks(k).x)-min(f_tracks(k).x))^2+(max(f_tracks(k).y)-min(f_tracks(k).y))^2);
        contourlength(1,k)=length(s_tracks(k).x);
        dev(1,k)=contourlength(1,k)/RG(1,k);
        if dev(1,k)<lengthdev
            ind_lengthdev(1,k)=1;
        end
    end
    trackstats.RG=RG;
    trackstats.contourlength=contourlength;
    trackstats.dev=dev;
else
    ind_lengthdev=ones(1,length(f_tracks));
    
end

if myoptions.minmaxcontourlength==1
    ind_minmaxcontourlength=zeros(1,length(f_tracks));
    if myoptions.lengthdev==0
        s_tracks=EqualSpace(f_tracks,1);contourlength=zeros(1,length(f_tracks));
        for k=1:length(s_tracks)
            contourlength(1,k)=length(s_tracks(k).x);
            if contourlength(1,k)> mincontourlength && contourlength(1,k)<maxcontourlength
                ind_minmaxcontourlength(1,k)=1;
            end
        end
        trackstats.contourlength=contourlength;
    else
        for k=1:length(contourlength)
            if contourlength(1,k)> mincontourlength && contourlength(1,k)< maxcontourlength
                ind_minmaxcontourlength(1,k)=1;
            end
        end
    end 
else
    ind_minmaxcontourlength=ones(1,length(f_tracks));
end

ind_total=ind_minmaxlengthtime & ind_minmedianspeed & ind_minRG & ind_lengthdev & ind_minmaxcontourlength;

f_tracks=f_tracks(ind_total);
end

function tracks_lambda=EqualSpace(tracks,spacing)
tracks_lambda=struct('x',[],'y',[],'s',[],'u',[],'v',[],'t',[]);
for i=1:length(tracks)
    s=zeros(1,length(tracks(i).x));
    for ii=2:length(tracks(i).x)
        s(ii)=s(ii-1)+sqrt((tracks(i).x(ii)-tracks(i).x(ii-1))^2+(tracks(i).y(ii)-tracks(i).y(ii-1))^2);
    end
    sq=0:spacing:s(end);
    x=interp1(s,tracks(i).x,sq);y=interp1(s,tracks(i).y,sq);u=interp1(s,tracks(i).u,sq,'pchip');v=interp1(s,tracks(i).v,sq,'pchip');t=interp1(s,tracks(i).t,sq);
    tracks_lambda(i,1).x(:,1)=x;tracks_lambda(i,1).y(:,1)=y;tracks_lambda(i,1).u(:,1)=u;tracks_lambda(i,1).v(:,1)=v;tracks_lambda(i,1).s(:,1)=s;tracks_lambda(i,1).t(:,1)=t;
end
end